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A model for crack propagation and contact line depinning is studied. Although the model con- 
tains nonlocal interactions, it obeys general scaling relations for depinning via localized bursts or 
avalanches. Our numerical result for the roughness exponent in one dimension, x ~ 0-49 ± 0.05, 
agrees with recent experiments on cracks measuring the in-plane roughness x — 0.5 — 0.6, as well as 
mean field arguments giving x ~ 1/2, but is significantly higher than the functional renormalization 
group prediction x = 1/3- 
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Despite its enormous practical importance and in- 
tense efforts to model the phenomenon, the physics of 
crack propagation is poorly understood [Q. Recently, 
Bouchaud et al. ^ proposed that crack propagation may 
be related to the dynamics of interface depinning by 
quenched random impurities Subsequently, Schmit- 
tbuhl Q et al. argued that the propagation of a crack in a 
certain idealized geometry has the same in-plane dynam- 
ics as contact line (CL) depinning on a dirty substrate. 
This mapping would apply specifically to a slowly ad- 
vancing crack confined within an easy plane between two 
elastic solids that are being slowly pulled apart. They 
numerically measured an in-plane roughness index x — 
0.35, in seeming agreement with the functional renormal- 
ization group (RG) 1^,^ prediction by Erta§ and Kardar 
§ that X = (2 - d)/3 exactly §. Da guier, Bouchaud, 
and Lapasset |^ have actually measured the in-plane 
roughness index of slowly moving cracks in metallic alloys 
and for stopped cracks found x — 0-5 — 0.6, in apparent 
contrast to the theoretical predictions. Since the model 
Q clearly represents a vastly simplified description of 
crack propagation, the discrepency might be due to an 
inherent deficiency. Here we present results which lead 
to a different conclusion. 

Our main numerical result from simulations of the 
model is X = 0.49 ± 0.05 with scaling over almost three 
decades in length scale. Our result is inconsistent with 
the functional RG prediction, but compares reasonably 
well with the experimental finding. We also show that de- 
spite the nonlocal interactions in the model motion near 
the depinning transition takes place in terms of localized 
bursts, or avalanches. As a result the critical dynamics 
obeys general scaling relations for avalanche phenomena 
in systems out of equilibrium |10|. Finally, based on nu- 
merical evidence, we conjecture that the critical expo- 
nents for CL depinning and the related model for crack 
propagation are given by a simple mean field theory giv- 
ing X = 1/2, D = 3/2, r = 4/3, = 2, 7 = 2, z = 1, 
df = 1/2, (3 — 1, and tt = 2 in one dimension. 



Applying a sufficient force to an interface causes it to 
move through a random medium with a finite velocity 
which vanishes continuously at a depinning transition 
where quenched impurities pin the interface into a static 
configuration. Rather than exhibiting a smooth continu- 
ous motion, the dynamics near the transition point takes 
place in terms of intermittent, localized bursts. A scal- 
ing theory has been developed that relates the critical 
dynamics of interface motion near the depinning transi- 
tion to the spatiotemporal structure of avalanches [p^ . 
This theory encompasses systems with local interactions 
[ ]lT| , p^ such as the motion of magnetic domain walls in 
the presence of quenched disorder, fluid invasion in a 
porous media, extinguishing flame fronts jist , invasion 
percolation, and flux creep |l^. Here we propose that 
it also includes contact line depinning and the related 
phenomena in crack propagation. 

Unlike the examples mentioned above, CL motion is 
governed by a nonlocal integral equation. This comes 
about because the capillary energy associated with small 
deformations of the CL has the long wavelength limit 
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associated with increasing the entire area of the liquid- 
vapor interface, rather than the usual line tension energy 
q2|/i(g)|2 ]i5|. Here h{x), the height of the CL profile 
at position x, is assumed to be single valued. O — >• 
is the macroscopic contact angle that the liquid-vapor 
interface makes with the substrate at the CL, 7 is the 
liquid-vapor interfacial tension, a is a small scale cutoff, 
and L is the linear extent of the one dimensional CL. 
The unusual capillary enery has a significant effect on 
the observed height fluctuations both in equilibrium and 
out of equilibrium []l5| . 

If the CL is driven slowly and most of the dissipation 
occurs in its vicinity, the capillary energy enters into the 
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equation of motion for the height profile as an applied 
force Fint{x,t) = — , where t is time. In addition, 
there is a random contribution to the local force density 
from the quenched impurities ryfx, /i(a;, t)) . The CL is 
driven with an external force F [M. 

In our numerical simulations the variables {h,x,t) are 
discretized to be integers. We impose a geometry where 
the CL is periodic in a;, i.e. h{x + mL) = h{x), where m 
is integer and L is the system size studied. This requires 
a summation of the nonlocal kernel over all periods and 
modifies the l/x^ interaction to be (^) / sin^ (ttx/L). In 
the thermodynamic L — > oo limit the summed interaction 
simplifies to l/x"^. The total force density f{x,t) 

— ) > „ „,x. ri[x,h) + F , (2) 
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where random forces rj are chosen independently from a 
flat distribution between zero and one, and the sum is 
over all lattice sites. Essentially the same equation was 
proposed |^ to describe the in-plane dynamics of a crack 
propagating within an easy plane of an elastic block that 
is being slowly wedged open. Of course, this represents 
a very idealized picture for crack propagation since, for 
instance, it ignores the out-of-plane meandering of the 
propagating crack and the out-of-plane roughness of the 
resulting crack surface 0]. 

A constant force depinning transition is implemented 
as follows: At each time step, the force at each site is 
evaluated. If f{x, t) > 0, then the site is unstable or 
"active" and the height at that site is advanced by one 
unit h — > h + I; otherwise the site is pinned and the 
height is unchanged. After the heights at all active sites 
have been advanced by one unit, time is advanced by one 
unit {t ^ t + 1), the local forces are re-evaluated, and 
the process is repeated. For F > Fc, the interface moves 
with a finite average velocity v = Uact/L where Uact is 
the average number of unstable, moving sites. 

The critical value Fc can be found by locating the site 
with the largest value of f{x, t) in a pinned configura- 
tion and increasing F until that site becomes unstable. 
This generates a burst of activity, or an avalanche, during 
which the interface moves. Eventually the burst dies out 
and the interface becomes stuck, with all sites frozen in 
a new configuration. The closer F gets to Fc, the larger 
is the average spurt of growth, or avalanche size, s. The 
average avalanche size diverges at the depinning transi- 
tion F = Fc- In order to obtain statistics for the steady 
state properties, one can set a value F — Fc — AF, where 
AF <C 1. At a moment in time when there are no ac- 
tive sites with f{x, t) > 0, the current F avalanche stops. 
Then a new F avalanche is initiated by advancing the 
height of the site with the largest value f{x,t) by one 
step. Starting from a flat configuration, h{x, t — 0) — 
for all X and repeating this process many times, the CL 
will reach a steady state in a time tss ~ L'^ ■ AH of the 



steady state properties of the depinning transition stud- 
ied in this work were obtained using this method. 

In Fig. 1, we show an actual realization of the height 
profile in the steady state near the depinning transition. 
The method of estimating the roughness exponent x 
which works best for this model, is to calculate the power 
spectrum P{k) for the height profile in the steady 
state. Fig. 2 shows power law behavior P{k) ~ k~^~^^. 
The asymptotic slope is —1.98, which fits the data over 
2.5 decades and gives the result x = 0.49. 

One can compare Fig. 2 to the equivalent Fig. 3. in 
Ref. Q where a slope for P{k) of —1.7 was measured with 
much fewer samples on a somewhat smaller system size 
L = 2048 than we study. On inspection of their figure, 
it is apparent that there is a systematic deviation from 
their fitted slope at small values of k (large values of x) 
where the asymptotic behavior dominates. Their data 
points at small k do not fall near the "best fit" line and 
actually give an apparent slope which is closer to what 
we measure. 

The equal time, height-height correlation function is 
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where the angular brackets denote an average over time t 
and over x, and h{t) is the average height of the interface 
at time t. In Fig. 3, we plot C{r)/L'^^ versus r/L, where 
X assumes its mean field value 1/2, so that C(G) ~ L. 
The large deviation for small system size L = 128 for r 
near zero indicates that such a small system size is not in 
the asymptotic regime. However we find reasonable good 
data collapse for L = 1024 — 8192. Our main numerical 
result is that x is bounded by x = 0.49 ± 0.05. 

We now discuss the scaling behavior of avalanches in 
the model. The size s of an F avalanche is the integrated 
motion, or the area between the initial configuration and 
the final configuration of the avalanche. In an experi- 
mental situation, perhaps, F would be increased slightly 
for each subsequent avalanche, but for a sufficiently large 
system, many avalanches would occur within a narrow in- 
terval centered on F. The statistics of those avalanches 
are described by the probability distribution P(s, AF). 

In analogy with other depinning problems |]l6| , ^ |lO | , it 
is plausible that the probability distribution of avalanche 
sizes is given by 



P(s,AF) - s-^g{sAF''^) 
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near the critical point AF = 0. Fig. 4 shows the mea- 
sured distribution for small AF; it decays as a power 
law over four decades with a characteristic exponent 
T — 1.31 ± 0.06 up to a cutoff determined by the sys- 
tem size. 

In spite of the non-local nature of the interactions, 
we find that similar to other interface models which 
have been studied 0,0,nfl, the dynamics during an 
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avalanche is local. Fig. 5 shows the location of the ac- 
tive, moving sites vs. time, for one avalanche of moderate 
size. It is clear from the figure that the avalanche spreads 
out in time. In fact, the distance spread r{t) grows with 
time as r{t) ~ t^^^. If the projection of the completed 
avalanche onto the sites that it covered is compact, and 
the avalanche size s ^ , where r is the spatial extent 
of the avalanche, then D — d + x In general D ^ z. 
For the same avalanches used in Fig. 4, we measured 
D = 1.5±0.15 which is consistent with the scaling rela- 
tion and our result for x- 

We now discuss a simple mean field theory to explain 
these results. The characteristic dimension of the non- 
local interaction term in Eq. (2) at length scale I is 
h/l where h ^ l^. The sum of random forces at this 
scale makes a contribution l^^^/l per unit length. Bal- 
ancing these two terms gives x = 1/2- Since the velocity 
v{l) ^ h/t f{l), the time i ^ Z so that the dynamical 
exponent z — 1. Defining the exponent f3 by v ^ AF^ 
and the cutoff Ico ^ AF~'^, one obtains the scaling rela- 
tion P = viz — x) The exponent /3 can be obtained by 
integrating the equation of motion over the entire system, 
so that the nonlocal part vanishes. This gives 



F+- dxr){x, h) ^ F + rj{F) 
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where the average ri{F) over the CL is evaluated at a 
particular value of F in the steady state. At F = Fc, 
V = so that 7y(Fc) = -F^. In general for F ^ Fc + AF 
we have ri{F) — —Fc + g{AF). Substituting this last 
result into Eq. (0), we get 



in the constant velocity and constant force case for this 
problem. 

In addition to CL depinning and crack propagation, 
exact results have been put forward using the functional 
RG for interface depinning with purely local interactions, 
where Narayan and Fisher predicted x = (4 — d) /3 [|| . 
Numerical simulations of related lattice models in both 
one and two dimensions give x = 1-2 3 |10[ (x = 1.25 
||l2|) in d = 1 and X = 0.72 in d = 2 [|10[. In all cases 
thus far where a comparison can be made, including our 
result, the numerical simulations give a roughness expo- 
nent higher than the functional RG prediction. 

Our theoretical picture for depinning phenomena is 
based on avalanche dynamics. To our knowledge, 
avalanches are not explicitly contained in the functional 
RG theory. We suspect that the apparent failure of the 
functional RG to reproduce numerical results for inter- 
face depinning may express an inadequacy of these meth- 
ods to describe avalanche dynamics. The agreement be- 
tween our numerical results and experiments measuring 
the in-plane roughness index of cracks supports the model 
for crack propagation proposed by Schniittbuhl, Roux, 
Vilotte, and Mal0y |||. 
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partment of Energy Division of Materials Science, under 
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partment of Energy Distinguished Postdoctoral Research 
Program for partial financial support. 
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Assuming g has a power series expansion gives (3=1 The 
scaling relation for /3 then gives i/ = 2 ||l^ . 

In general, for systems with local avalanche dynam- 
ics the exponent r for the distribution of avalanche sizes 
obeys the relation r = and substituting the 

mean field values for one dimension {d = 1) gives r = 4/3 
in agreement with our numerics. The exponent 7 for the 
divergence of the average avalanche size is 7 = i^D(2 — r), 
and 7 = 2 in mean field. In the constant velocity depin- 
ning transition the critical exponent b in Ref. [Q charac- 
terizing the power law distribution between subsequent 
activity obeys b = tt — 1 = D{2 — t) according to the 



scaling theory ]13|,|lC| 
be 6 = 0.9 ± 0.05 
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FIG. 1. An actual realization of the height profile h{x,t) 
for L = 8192 [20]. 

FIG. 2. A log-log plot of the power spectrum of the height 
profile versus k for L = 8192. The height profile was sampled 
every t = 273 time steps and over 10^ samples were used to 
compute P{k). The straight line is a best fit to all the points 
k < 0.1. It has a slope of -1.98. 



FIG. 3. Rescaled plot of the equal time height-height cor- 
relation function (C(r)/L) vs. r/L for L = 128, F = 0.83 
(open circles), L = 1024, F = 0.828 (triangles), and L = 8192 
(filled circles). This collapse of the data for L = 1024 - 8192 
is consistent with x = 1/2. 

FIG. 4. The avalanche size distribution near the depinning 
transition for L = 8192. Over 5 x 10* avalanches were mea- 
sured. The slope gives an estimate of t = 1.31 ± 0.06. 

FIG. 5. The location (marked as dots) of the active, mov- 
ing sites as a function of time, for L = 1024, during an an 
avalanche which begins at t = and ends at t = 545. The 
entire avalanche has a size s = 9814. The avalanches are local 
and spread out in time as r ~ t^^^ . 
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